﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "gjfq")]
    public static unsafe void gjfq(double t, double h, int n, IntPtr y_ptr, double eps, int k, IntPtr z_ptr, IntPtr f_x_ya_n_da_ptr)
    {
        double* y = (double*)y_ptr.ToPointer();
        double* z = (double*)z_ptr.ToPointer();
        f_x_ya_n_da = Marshal.GetDelegateForFunctionPointer<delegatefunc_x_ya_n_da>(f_x_ya_n_da_ptr);

        gjfq(t, h, n, y, eps, k, z);
    }

    /// <summary>
    /// 全区间积分双边法
    /// f计算微分方程组中各方程右端函数值的函数名。
    /// void gjfq(double t, double h, int n, double y[], double eps, int k, double z[], void (*f)(double,double [],int,double []))
    /// </summary>
    /// <param name="t">积分起始点。</param>
    /// <param name="h">积分步长。</param>
    /// <param name="n">一阶微分方程组中方程个数，也是未知函数个数。</param>
    /// <param name="y">y[n]存放n个未知函数在起始点t处的函数值。</param>
    /// <param name="eps">变步长Runge_Kutta法的控制精度要求。</param>
    /// <param name="k">积分步数（包括起始点这一步）。</param>
    /// <param name="z">z[n][k]返回k个积分点（包括起始点）上的未知函数值。</param>
    public static unsafe void gjfq(double t, double h, int n, double* y, double eps, int k, double* z)
    {
        int i, j;
        double a, qq;
        double* d = stackalloc double[n];
        double* p = stackalloc double[n];
        double* u = stackalloc double[n];
        double* v = stackalloc double[n];
        double* w = stackalloc double[n];

        for (i = 0; i <= n - 1; i++)
        {
            p[i] = 0.0;
            z[i * k] = y[i];
        }
        a = t;
        f_x_ya_n_da(t, y, n, d);
        for (j = 0; j <= n - 1; j++)
        {
            u[j] = d[j];
        }
        runge_kutta(t, h, n, y, eps);
        t = a + h;
        f_x_ya_n_da(t, y, n, d);
        for (j = 0; j <= n - 1; j++)
        {
            z[j * k + 1] = y[j];
            v[j] = d[j];
        }
        for (j = 0; j <= n - 1; j++)
        {
            p[j] = -4.0 * z[j * k + 1] + 5.0 * z[j * k] + 2.0 * h * (2.0 * v[j] + u[j]);
            y[j] = p[j];
        }
        t = a + 2.0 * h;
        f_x_ya_n_da(t, y, n, d);
        for (j = 0; j <= n - 1; j++)
        {
            qq = 2.0 * h * (d[j] - 2.0 * v[j] - 2.0 * u[j]) / 3.0;
            qq = qq + 4.0 * z[j * k + 1] - 3.0 * z[j * k];
            z[j * k + 2] = (p[j] + qq) / 2.0;
            y[j] = z[j * k + 2];
        }
        for (i = 3; i <= k - 1; i++)
        {
            t = a + (i - 1) * h;
            f_x_ya_n_da(t, y, n, d);
            for (j = 0; j <= n - 1; j++)
            {
                u[j] = v[j];
                v[j] = d[j];
            }
            for (j = 0; j <= n - 1; j++)
            {
                qq = -4.0 * z[j * k + i - 1] + 5.0 * z[j * k + i - 2];
                p[j] = qq + 2.0 * h * (2.0 * v[j] + u[j]);
                y[j] = p[j];
            }
            t = t + h;
            f_x_ya_n_da(t, y, n, d);
            for (j = 0; j <= n - 1; j++)
            {
                qq = 2.0 * h * (d[j] - 2.0 * v[j] - 2.0 * u[j]) / 3.0;
                qq = qq + 4.0 * z[j * k + i - 1] - 3.0 * z[j * k + i - 2];
                y[j] = (p[j] + qq) / 2.0;
                z[j * k + i] = y[j];
            }
        }
        return;
    }
    /*
    // 全区间积分双边法例
      int main()
      { 
          int i,j;
          void  gjfqf(double,double [],int,double []);
          double y[2],z[2][11];
          double t,h,eps;
          t=0.0; h=0.1; eps=0.0000001;
          y[0]=1.0; y[1]=0.0;
          gjfq(t,h,2,y,eps,11,&z[0][0],gjfqf);
          for (i=0; i<=10; i++)
          { 
              t=i*h;
              cout <<"t = " <<t;
              for (j=0; j<=1; j++)
                  cout <<"  y(" <<j <<") = " <<setw(10) <<z[j][i];
              cout <<endl;
          }
          return 0;
      }
    // 计算微分方程组中各方程右端函数值
      void gjfqf(double t, double y[], int n, double d[])
      { 
          t=t; n=n;
          d[0]=-y[1]; d[1]=y[0];
          return;
      }
    */
}

